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Abstract. Following on after three previous papers discussing the formation of 
primordial black holes during the radiation-dominated era of the early universe, 
we present here a further investigation of the critical nature of the collapse. In 
particular, we focus on the long-lived intermediate state, which appears in collapses of 
perturbations close to the critical limit, and examine the extent to which this follows 
a similarity solution, as seen for critical collapse under more idealized circumstances 
(rather than within the context of an expanding universe, as studied here). We find that 
a similarity solution is indeed realised, to good approximation, for a region contained 
within the past light-cone of the forming black hole (and eventual singularity). The 
self-similarity is not exact, however, and this is explained by the presence within the 
light-cone of some outer matter still coupled to the expanding universe, which does not 
participate in the self-similarity. Our main interest, from a cosmological point of view, 
is in a radiative fluid with equation of state parameter w = 1/3 (when the pressure 
p and energy density e are taken to be related by p — we). Other values of w, in 
the range — 1, have also been considered in the literature on critical collapse and we 
have looked at some of these too, within the context of our calculations, with the aim 
of gaining further insight into our main case of interest. As expected, we find that 
the features of scaling-law behaviour, intermediate state and similarity solution arc 
preserved in all of the cases studied but with some interesting variations in the details. 
As in our previous work, we have started our simulations with initial supra-horizon 
scale perturbations of a type which could have come from inflation. 
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1. Introduction 

A population of primordial black holes (PBHs) might have been formed in the early 
universe due to gravitational collapse of large-enough cosmological perturbations. 
Following the initial papers suggesting this (Zel'dovich & Novikov (1969) [JJ; Hawking 
(1971) [2]; Carr [31 H]), various authors then investigated the process numerically 
(Nadezhin, Novikov & Polnarev (1978) [5]; Bicknell & Henriksen (1979) ©; Novikov 
& Polnarev (1980) [?J Niemeyer and Jedamzik (1999) [8]; Shibata and Sasaki (1999) [9]; 
Hawke & Stewart (2002) [TO]; Musco, Miller & Rezzolla (2005) p]). PBHs would have 
been formed if some cosmological perturbations had reached an amplitude 5, at the time 
of cosmological horizon crossing, greater than a certain threshold value S c (the S is often 
defined as the relative mass excess inside the overdense region, measured at the time of 
horizon crossing). The simulations clarified many different aspects of the nature of the 
collapse, with particular reference to the radiative era of the universe (with the equation 
of state being taken as p = |e, where p is the pressure and e is the energy density). 
This context represents the most commonly studied scenario for PBH formation. 

Niemeyer & Jedamzik (1999) [HI H2] showed that the masses of PBHs produced in 
the radiative era by perturbations of a given profile type, follow the typical scaling-law 
behaviour of critical collapse, first discovered for idealized conditions by Choptuik (1993) 
[Ej, i.e. the masses of the black holes produced follow a power law Mbh oc (5 — 5 C ) 7 
if 5 is close enough to 5 C . Neilsen & Choptuik (2000) [15] later showed that one could 
get a critical collapse, using a succession of "imploding shells" of matter as the initial 
conditions, for p = we equations of state with any w in the range — 1. For these 
configurations, they found that the value of the critical exponent 7 was dependent only 
on the value of w and not on the particular form of the perturbation profile. 

In 2002, Hawke & Stewart [10] came back to the problem of critical collapse in the 
context of PBH formation in the early universe with a radiative fluid, and investigated 
the nature of the collapse going down to smaller values of (5 — 5 C ) than Niemeyer 
& Jedamzik had been able to do with their code. For the larger values of (5 — 5 C ), 
comparable to those of [8], Hawke & Stewart again found a scaling law with a similar 
value of 7, but for smaller (5 — 5 C ) they saw formation of strong shocks and the curve 
flattened off at a minimum mass of around 10~ 3 of their horizon mass. They also found 
that the value of S c depended strongly on the shape of the perturbation, in contrast to 

i- 

In Musco et al (2005) [UJ, we considered PBH formation with initial conditions 
given by linear perturbations of the energy density and velocity fields approximating 
the growing components of cosmological perturbations, and imposed with a length- 
scale much larger than the cosmological horizon. The domination of the growing 
component becomes even greater by the time of horizon crossing, as any residual 
decaying component dies away. We again saw a scaling law with similar 7 for the 
range of values of (5 — 5 C ) used by [5] but found very different values of 5 C from [S] 
for similar perturbation shapes (we were concentrating on those types of perturbation). 
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This difference was attributed to the fact that in [8J, the perturbations were made just 
in the energy density and were imposed directly at the horizon scale, so that their value 
of S c was calculated with inclusion of a substantial decaying component which did not 
then contribute to the black hole formation. If one focuses on the effect of perturbations 
originating from inflation, then it is clearly only growing components which are relevant 
at much later times. 

In Polnarev & Musco (2007) [16J, this kind of cosmologically-relevant initial 
perturbation was imposed in a more precise way, using an asymptotic quasi- 
homogeneous solution [17J. Starting from a curvature perturbation, which is a time- 
independent quantity when the perturbation length-scale is much larger than the 
cosmological horizon [IS] , perturbations in all of the other quantities can then be 
specified in a consistent way, giving a solution with only a growing component. The 
metric perturbation can be large even when the perturbations in energy density and 
velocity are small, and only when there is a large-amplitude metric perturbation, 
corresponding to a non-linear initial perturbation of the curvature, can one get 5 > 5 C 
at horizon crossing. 

In Musco et al (2009) [IS], we returned to the issue raised by Hawke & Stewart [TU] 
concerning whether the scaling law would continue down to very small values of (5 — 5 c ). 
In order to address this, we needed to modify our code with the inclusion of adaptive 
mesh refinement (AMR) so as to be able to handle the extreme conditions which arise 
near to the critical limit. In view of our earlier work, we were particularly wanting 
to investigate the effect of using initial perturbations with just a growing component 
imposed outside the horizon scale, rather than the non-linear sub-horizon scale initial 
perturbations used in JTU] . We again used the quasi-homogeneous solution to provide our 
initial conditions. Doing this, we did not observe the shock formation seen in [10] (which 
we attributed to the presence of a non-linear decaying component in their calculations) 
but instead got a regular scaling law behaviour all the way down to the vicinity of 
the resolution limit of our scheme (going beyond the most extreme values shown in 
[10J). A striking feature of our calculations was the appearance of an "intermediate 
state" during collapse for cases near to the critical limit: the collapse proceeded to a 
compactness 2M/R ~ 0.5 (where R and M are the current radius and mass of the 
condensation) which then remained roughly constant for a large number of dynamical 
timescales, with the condensation progressively shedding matter via a relativistic wind 
and shrinking, until eventually it either reached a mass at which the final collapse to a 
black hole could take place or started to disperse. 

The existence of a similarity solution, for collapses near to the critical limit, is a 
key feature of the theory of critical collapse [201 EI], arid in [19] we saw some evidence of 
self-similarity also in the context of PBH formation within an expanding universe. This 
is associated with the long-lived intermediate state. In the present paper, we return to 
investigate this further, focusing on the intermediate state and examining the extent to 
which it does, in fact, follow a similarity solution. We also investigate other values of w. 
While our main interest from a cosmological point of view is in w = 1/3, it is useful to 
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look at other values as well, using the same approach, so as to get further insight into 
the main case of interest. Also, other values of w might be relevant for PBH formation 
at certain particular times in the early universe. 

The results presented in this paper have been obtained using the same numerical 
code as in [TJJ] , with some fine tuning of the AMR and modifications for calculating the 
past light cone of the forming black hole. After this Introduction, Section 2 reviews our 
mathematical formulation of the problem, Section 3 contains a brief summary of the 
numerical methods used and also a description of our newly-introduced calculation of 
the past light cone; Section 4 presents results from our investigations of the intermediate 
state and the similarity solution, and of the consequences of using different values for 
w; Section 5 gives conclusions. In an Appendix, we discuss the concepts of cosmological 
and black-hole horizons within the present treatment. Throughout, we use units for 
which c = G = 1, except where otherwise stated. 

2. Mathematical formulation of the problem 

2.1. Cosmic-time slicing 

For the calculations described here, we have followed the same basic methodology as 
described in our previous papers [HJ [161 HH] (which we will refer to as Papers 1, 2 and 
3 respectively). We therefore give just a brief summary of it here; more details are 
contained in the previous papers. 

We use two different formulations of the general relativistic hydrodynamic 
equations: one for setting the initial conditions and the other for studying the black 
hole formation. Throughout, we are assuming spherical symmetry and that the medium 
can be treated as a perfect fluid; we use a Lagrangian formulation of the equations with 
a radial coordinate r which is co-moving with the matter. 

For setting the initial conditions, it is convenient to use a diagonal form of the 
metric, with the time coordinate t reducing to the standard Friedmann-Robertson- 
Walker (FRW) time in the case of a homogeneous medium with no perturbations. (This 
sort of time coordinate is therefore often referred to as "cosmic time"). We write this 
metric in the form given by Misner Sz Sharp (MS) (22], whose approach we follow here 
in writing the GR hydrodynamic equations: 



with the coefficients a, b and R being functions of r and t and with R playing the role 
of an Eulerian radial coordinate. Using the definitions 



ds 2 = -a 2 dt 2 + b 2 dr 2 + R 2 (d6 2 + sin 2 9dcp 2 ) 



(1) 




(2) 




(3) 
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(4) 



one defines the quantities 

U = D t R , 
and 

r = D r R , (5) 

where U is the radial component of four-velocity in the "Eulerian" frame and T is a 
generalized Lorentz factor. The metric coefficient b can then be written as 

With these specifications, the GR hydrodynamic equations can be written in the 
following form (with the notation that e is the energy density, p is the pressure, p 
is the compression factor and M is the mass contained inside radius R): 



D t U- 

D t p = 

D t e = 

D t M 
D r a = 



r m 

_{e + p) R 2 



P 



TR 2 

e + p 



D r (R 2 U), , 



P 



Dtp, 



-4ttR 2 pU , 
a n 

—D r Pi 

e + p 



D r M = 47iR 2 Te , 
plus a constraint equation 



r 2 = i + u 2 



2M 



(7) 
(8) 

(9) 
(10) 
(11) 
(12) 

(13) 
(14) 



An equation of state is also needed, and we are here considering ones of the form 
p = we , 

where w is a constant parameter which can take various values. When a cosmological 
constant is added, the total energy density e and the total pressure p are given by 

A ~ A 

e = e+— p = p -—- (15) 

671 07T 

where A/87T is the additional vacuum component of the energy density and —A/8tt is 
the vacuum component of the pressure; these vacuum contributions satisfy (14) with 
w — — 1. The cosmological constant produces a gravitational effect that changes the 
background, but does not affect the fluid in other ways. 
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2.2. Initial conditions 

The initial conditions are set by introducing a perturbation of the otherwise uniform 
medium representing the cosmological background solution, with the length-scale of the 
perturbation R being much larger than the cosmological horizon Rh = H -1 . Under 
these circumstances, the perturbations in e and U can be extremely small while still 
giving a large-amplitude perturbation of the metric (as is necessary if a black hole 
is eventually to be formed) and the above system of equations can then be solved 
analytically to first order in the small parameter e = (Rh/Ro) 2 < 1- A full discussion 
of this has been given in Paper 2. The result obtained, referred to as the "quasi- 
homogeneous solution", gives formulae for the perturbations of all of the metric and 
hydrodynamical quantities in terms only of a curvature perturbation profile K(r) that 
is conveniently time-independent when e C 1 [18J. 

To characterize the amplitude of the perturbation, we use the integrated quantity 

which measures the relative mass excess within the overdense region, as frequently done 
in the literature; e& is the background value of the energy density and ro is the co-moving 
length-scale of the overdense region of the perturbation, determined by the particular 
expression used for K(r). In Paper 2 (to which we refer for details) it was shown that 
e(r,t)-e»(t) 3(1 + w) rg d [r 3 K(r)] 

e b (t) U 5 + 3w 3r 2 dr { ' 

and here, as in Paper 3, we make the simplest choice of using a Gaussian for the initial 
curvature profile K(r), 

K{r)=exp(--^\ , (18) 

which implies rg = 3A 2 and leads to the corresponding Mexican Hat profile for the 
energy density at the initial time: 

(19) 

(.xpl-^(-) I, (20) 




3(1 + 10) A2 



1 - 



r 
Jo, 

where A is the free parameter used at the initial time to vary the perturbation amplitude. 



= e(t) v 7 A 
e b (t) W 5 + 3w 




Inserting (20) into (Il6[) one gets 



9(1 + 10) . 2 ( 3 



^ = ^TW Aex H"2j' (21) 

within the linear regime where e C 1, with the time evolution of the perturbation being 
given only by e(t), which measures how the cosmological horizon Rjj grows with respect 
to the length-scale of the perturbation Rq. Inserting the expression for e(t) gives 

2(l+3iu) 

m = mU) , (22) 
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which leads to the familiar relationships 5(t) oc t when w — 1/3 and 8(t) oc t 2 ^ 3 when 
w — 0. Full details about setting initial conditions using the quasi-homogeneous solution 
are given in Paper 2. 

For our subsequent discussion of the scaling laws, we need a value for 5 measured 
at a standard time and it is convenient to take this as being the time at which e(t) = 1. 
For small perturbations, this coincides with the horizon-crossing time, and it does not 
differ by very much from that even for the larger perturbations of interest for PBH 
formation, at least for the type of perturbation profile being used here. Since evaluating 
(21) analytically is much more precise than making a numerical integral on the grid at 
horizon-crossing time, we use that value for the discussion of the scaling laws, as we did 
in Papers 2 and 3. 

Profiles similar to a Mexican Hat for the energy density perturbation represent a 
plausible and simple choice, and previous studies [13l HU HE] have shown that other 
centred peak perturbation profiles, such as polynomial ones, behave in a very similar 
way, giving almost the same values for 5 C . Off-centred profiles could, in principle, give 
rise to more substantial changes but are rather unlikely to arise naturally. Therefore, in 
the absence of any strong motivation for considering those kinds of profile, we restrict 
our choice of initial conditions only to Mexican Hat profiles given by (20). 



2.3. Observer-time slicing 

The MS approach using cosmic time slicing is convenient for setting initial conditions but 
has the well-known drawback for calculations of black hole formation that singularities 
are formed rather quickly when using it and further, potentially observable, evolution 
cannot then be followed unless an excision procedure is used. Various slicing conditions 
can be used to avoid this difficulty but for calculations in spherical symmetry it is 
particularly convenient to use null slicing. In our work we have used the "observer 
time" null-slicing formulation of Hernandez & Misner [2l] where the time coordinate is 
taken as the time at which an outgoing radial light ray emanating from an event reaches 
a distant observer. (In the original formulation, this observer was placed at future null 
infinity but for calculations in an expanding cosmological background we use an FRW 
fundamental observer sufficiently far from the perturbed region so as to be unaffected by 
the perturbation.) We use the MS approach for setting up initial data and for evolving 
it so as to produce data on a null slice. This is then used as input for our observer-time 
code with which we follow the black hole formation. 

For the observer-time calculation, the metric ([I]) is re- written as 

ds 2 = -f 2 du 2 - 2fb dr du + R 2 (d6 2 + sin 2 6d<p 2 ) , (23) 

where u is the observer time and / is the new lapse function. The operators equivalent 
to @ and Q are now 

1 / d \ 
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(25) 



where D k is the radial derivative in the null slice and the corresponding derivative in 
the Misner-Sharp space-like slice is given by 

D r = D k - A • (26) 

The hydrodynamic equations can then be formulated in a way analogous to what was 
done in cosmic time. The observer-time equations used in the code, replacing the 
cosmic-time ones Q - (12), are: 

1 



D+U 



r m 

i 2 (~^-^ Dk V + ™ + ^Rp 
1 - c% [ (e + p) R 2 



+ c 2 s \D k U + 



2UT \ 
~1T 



Dtp 



Ae 



D t U - D k U 



2UT 



e + p 
P 

F + U) 
f 



-4irR(e + p)f . 



D k M = A%R 2 [eT — pU] , 



(27) 
(28) 

(29) 

(30) 
(31) 



where c s = a/ (dp/de) = y/w is the sound speed. The quantity T is given by (13), as 
before, and we also have 

T = D k R-U, (32) 

which replaces ^ in this slicing. 

3. The method of calculation 

The present calculations have been made with the same code as used in Paper 3. Since 
this has been fully described previously, we will just give a brief outline of it here. It 
is an explicit Lagrangian hydrodynamics code based on that of Miller & Motta (1989) 
[25] but with the grid organized in a way similar to that of Miller & Rezzolla (1995) 
[26] which was designed for calculations in an expanding cosmological background. The 
code has a long history and has been carefully tested in its various forms. The basic 
grid uses logarithmic spacing in a mass-type comoving coordinate, allowing it to reach 
out to very large radii while giving finer resolution at small radii. Our initial data is 
derived from the quasi-homogeneous solution and is specified on a space-like slice (at 
constant cosmic time) with e = 10 -2 , giving i? = 10 Rh- The outer edge of the grid 
has usually been placed at 90 Rh, which is far enough away so that there is no causal 
contact between it and the perturbed region during the time of the calculations. The 
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initial data is then evolved using the MS equations ((7 13), so as to generate a second set 
of initial data on a null slice (at constant observer time) and the null-slice initial data 
is then evolved using the observer-time equations [TT] . 

For the calculations presented in Paper 3, we introduced an adaptive mesh 
refinement scheme (AMR), on top of the existing logarithmic grid, giving us sufficient 
resolution so as to be able to follow black hole formation down to extremely small values 
of (5 — 5 C ). Having the AMR is particularly important for allowing us to follow the deep 
voids which form outside the central collapsing region in cases very close to the critical 
limit. The same scheme has been used with just minor modifications for the calculations 
presented here. Our aim in writing the AMR was to avoid the use of artificial viscosity, 
but it has now emerged that some residual artificial viscosity was still present and is 
necessary for correct functioning of the code in its present form. The presence of this 
is not thought to affect the results presented in any important way, however. We have 
successfully used the scheme with more than thirty levels of refinement and all relevant 
features of the solutions have been fully resolved. 

For the discussion of the similarity solutions in the next section, we are needing to 
know the location of the past light cone of the forming black hole (and of the eventual 
singularity). This involves tracing the paths of ingoing radial light rays. Considering 



the observer-time metric (23): along radial null rays, we have ds = d9 = dip = 0, so 
that 



fdu 2 + 2fbdrdu = 0, (33) 



giving 



fdu(fdu + 2bdr) = . (34) 

This has two solutions: du = 0, referring to an outgoing radial null ray (the basic 
definition for the observer-time as the time at which events are seen by a distant 
observer), and dr = — (f/2b)du for an ingoing null ray. The general expression for 
changes in R along a radial worldline is 
dR OR 

dR = —du + —dr . (35) 
ou or 

Using Q, Q and Q - (@ one has 

~^ = fU and ^=b(T + U), (36) 



and inserting these into (35) gives 

dR = fUdu + b(T + U)dr } (37) 
which, with the relation dr = —(f/2b)du, leads to 

f u = {(U-T), (38) 

along an ingoing radial null ray. 

The use of observer time gives a particularly clear treatment of light cones. We are 
interested in the past light cone of the asymptotically-forming black hole (which would 
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also then be the lightcone of the subsequent singularity). Our code does not allow us 
to make an evolution backwards in time, however, and so we need to find the location 
of this light cone at our initial time by means of an iterative process, looking for the 
ingoing light-ray which settles into the surface of the black hole at formation. We used 
a secant iteration for doing this. 

4. Results: intermediate state, similarity solutions and variation of w 

In this section, we present results from our further investigation of the intermediate state 
including, in particular, examining the extent to which it follows a similarity solution. 
We also investigate the effects of using values of w different from 1/3. As mentioned 
previously, while our main interest from a cosmo logical point of view is in w — 1/3, it 
is useful to look at other values as well, using the same approach, so as to gain further 
insight into the main case of interest. Some of these other values might also be directly 
relevant at some particular stages in the early universe. 

After a growing perturbation enters within the cosmological horizon, it will start to 
contract and then collapse if its amplitude 5 is large enough. If 5 is larger than the critical 
value 8 C , it will go on to produce a black hole, while for smaller values it will eventually 
disperse into the background medium. For 5 close to 5 C (either above or below it) one 
finds that the compactness 2M/R of the collapsing material first decreases, but then 
attains a value at which it remains for many dynamical timescales. This is the long-lived 
"intermediate state" and, for a radiation fluid, this has 2M/R ~ 0.5, quite far from the 
condition for a black hole. In this state, it sheds matter via a relativistic wind, which 
is driven by a huge pressure gradient at its surface, and a deep void opens up between 
it and the rest of the universe. As it loses matter via the wind, it shrinks maintaining 
an almost constant value of 2M/R. Eventually, its mass decreases to a point where it 
begins to deviate from the intermediate state, either collapsing to form a black hole (if 
5 > 8 C ) or dispersing (if 5 < S c ). This process was investigated in some detail in Paper 3 
for a radiation fluid with w — 1/3 and some evidence was seen of self-similar behaviour 
during the period of the intermediate state. Here, we return to this with the aim of 
answering two main questions: (i) Does the intermediate state really follow a similarity 
solution and how does this fit with the surrounding expanding universe? (ii) What is the 
effect of varying w, and do we find results for this in accordance with earlier calculations, 
where critical collapse was being studied under simple idealized circumstances rather 
than following cosmological perturbations within an expanding universe? 

When calculations are made, using our approach, for non-standard values of w, 
there is some artificiality because of having the whole universe following an unusual 
equation of state and also imposing initial conditions under those circumstances. If 
one were really considering in detail possible PBH formation with a variant equation 
of state, one would need a more sophisticated set-up than our present one. However, 
our aim here is just to get some general indications about the effect of varying w and 
to make a link with the previous work on critical collapse. We have studied a range of 
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Figure 1. In the top panels the behaviour of 2M/ R is plotted against R/Rh for 
a nearly-critical case (5 ~ 5 C ) for w = 0.1 and w = 1/3, at different time levels. 
Rh is the cosmological horizon scale at the moment of horizon crossing. The 
dashed curve indicates the initial conditions used by the observer-time code. 
The bottom panels show the time evolution of the peak of 2M/R during the 
intermediate state: tjtu is the observer time measured in units of the horizon 
crossing time, while (A) and (B) represent the limits of the time interval used 
in figure [2j dealing with the self similar solution. 



values for w between 0.01 and 0.6, looking at the intermediate state and the scaling law 
in each case. As w is increased, the features of the relativistic wind and the opening 
up of the void become progressively more extreme, and this eventually becomes very 
challenging for the adaptive scheme. It is because of this that we have just used values 
of w < 0.6 for the simulations presented here. 

We first compare the intermediate state for the standard case of w — 1/3 with that 
for a representative lower value, w = 0.1 (and then comment later on what happens for 
other values). The time spent in the intermediate state becomes progressively longer 
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as 5 approaches 5 C . Figure [T] shows the behaviour of 2M/R plotted against R/Rh 
at different time levels when a black hole is formed with a value of 5 extremely close 
to the threshold (Rh being the cosmological horizon scale at the moment of horizon 
crossing). The dashed curve shows the initial conditions used by the observer-time 
code, after which the evolution proceeds in the direction of decreasing {2M / R) peak until 
the intermediate state is reached, after which the peak moves progressively inwards. 
When the collapse is extremely close to the critical limit, as it is here, the intermediate 
state lasts almost up until the formation of the black hole, when (2M / 'R) pe ak departs 
from its almost constant value in the intermediate state and increases up to the black- 
hole value of 1. Another view of this behaviour is shown in the bottom panels of figure 
[TJ where the time evolution of (2M/ R) pe ak is plotted, for the same two cases as shown 
in the top panels, starting when the collapse is getting into the intermediate state. One 
can see that, for the smaller value of w, the intermediate state lasts for longer, with 
(2M / 'R) P eak being more constant. Note that in both cases the final collapse to form the 
black hole, when it comes, is extremely rapid in comparison with what has gone before 
(this is a characteristic feature when 5 is close to 5 C ). The behaviour for the other values 
of w which we have tested follows the trends seen for the two cases shown above. We 
will comment more about this in connection with figure [4] below. 

Turning now to the self-similarity, in figure [2] we plot quantities similar to ones 
used by Evans & Coleman [20] in their study of critical collapse and self-similarity: 
2M/R, AnR 2 e and U/c are plotted against £ = [R/C(t c — t) n ] at a succession of output 
times between the points marked with (A) and (B) in figure [2j Here, n is a self-similar 
exponent, C is a constant (= Raft™) which makes £ dimensionless, and t c is the critical 
time at which the singularity would be formed in the limiting case 5 = 5 C (we find 
t c in the way described in Paper 3). The time ordering of the curves is from top to 
bottom on the right hand side for 2M/R and 4nR 2 e, and from bottom to top in the 
mid section for U/c. Evans & Coleman [20] found n ~ 1.15 for the circumstances which 
they were studying (for a radiation fluid but with a different time-slicing from ours). 
We calculated best-fit values of n for our calculations and found n = 0.90 for w = 1/3 
and n = 0.83 for w = 0.1. The self-similarity seen in figure [2] is clearly not exact but 
represents a rather good approximation for all of the collapsing region out to the peak of 
2M/R and somewhat beyond. (Note that in the frames for w = 1/3, we have purposely 
included a final time-level where the self-similarity is just starting to break.) Further 
out, the solutions undergo a transition into the region of the wind and the void where a 
connection still remains with the surrounding universe. The influence of the connection 
between the collapsing matter and the surrounding universe is made clear by considering 
the past light cone of the forming black hole and the eventual singularity (calculated as 
described in Section 3). The location of the light cone is shown by the vertical dashed 
lines in figure [2] (there is a slight inward motion of this with respect to £ during the time 
shown, indicated by the gap between the two dashed lines). The self-similar region is 
completely contained within the light cone which, however, also contains material not 
participating in the self-similarity which can then stop the self-similarity being exact. 
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Figure 2. Self-similar behaviour during the intermediate state. For the 
two cases shown in figure [TJ we plot 2M/R, AirR 2 e and U/c against £ = 
[R/C(t c — t) n ] at a succession of output times between the points marked with 
(A) and (B) in figure [TJ The best-fit values for the exponent n are 0.90 for 
w = 1/3 and 0.83 for w = 0.1. The vertical dashed lines delimit the small 
range of variation of the location of the past light-cone (in terms of £) during 
this time; it is slightly moving inward. 
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Figure 3. Scaling behaviour for Mbh as function of (5 — S c ) for values of w 
ranging from 0.1 to 0.6; Mh is the cosmological horizon scale at the moment 
of horizon crossing. The best-fit values for 7 are indicated in each case. 
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The last outputs included in figure [2] are ones just before the deviations away from 
the approximate self-similarity in the inner regions start to become significant. In the 
left-hand column, we have actually included one last output where the deviations do 
start to be clearly seen, as an indication of the behaviour. Note that the breaking of 
the similarity occurs a little before the final dynamical collapse commences. 

In Figure [3] we show the scaling- law behaviour of Mbh plotted as a function of 
(5 — 5 C ) for several values of w ranging from 0.1 to 0.6. Good scaling laws are obtained 
in each case over the range of values of (5 — 8 C ) that we were able to study (note 
that this range was much smaller in the cases with w > 1/3 because of the extreme 
conditions encountered there). The best-fit values of 7 are completely consistent with 
those obtained previously by Maison [27] and by Neilsen and Choptuik [15] with an 
imploding shells calculation. 

In the left hand panel of figure [4], we show results obtained with our code for the 
dependence of 5 C on w. Very similar behaviour to this is obtained also for {2M / R) peak 
during the period of the intermediate state, but with (2M/R) pea k being slightly larger 
than 8 C . The first estimation of the value of 5 C was made in the pioneering paper by 
Carr [3] in 1975, using an approximate argument based on the Jeans length. He found 
8 C ~ w and our numerical simulations are in reasonable agreement with that, bearing in 
mind the approximations. This relation between 8 C and w confirms that any epochs in 
the early universe when the equation of state softens would be favorable for enhancing 
PBH production. One example of this could be the QCD phase transition, which has 
previously been discussed in connection with a PBH model for MACHOs [13J. In the 
right hand panel of figure [| we plot our values for 7 against w, for values of w up to 0.6 
(marked with crosses), and this shows a regular behaviour consistent with the limit of 




Figure 4. The left-hand plot shows the dependence of 5 C on the value of w, 
and the right-hand plot shows the behaviour of 7 as a function of w (which is 
almost linear within the range studied here). 
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7 — > 0.106 for w obtained by Snajdr [28]. The values of 7 found in our calculations 
are almost indistinguishable on this scale from those obtained in J27J EE] ■ 

Finally we want to comment on one of our earlier results, in Paper 1, where we 
investigated critical collapse for a radiative fluid also in the presence of a cosmological 
constant A, as described in section 2 here by equations (14) and (15). At the time of 
Paper 1, the code was not able to get very close to the critical limit because of not 
yet having the AMR, and the lowest value of (5 — 5 C ) that we were able to treat was 
around 10~ 3 . In that regime we found that the presence of a cosmological constant was 
affecting the scaling law, giving a change in the value of 7 in the sense of decreasing it 
with increasing (positive) A. Now, with the present AMR code, we have made similar 
calculations going down to much smaller values of (5 — 5 C ) and have seen that this 
change in 7 disappears as one gets closer to 5 C (i.e. the gradients of the two scaling 
laws converge to the same value). The reason for this is clear: if 5 is very close to 
5 C , the mass of the black hole is small and the fluid densities involved in forming it 
are high. Under these circumstances, the energy density related to the cosmological 
constant becomes negligible in comparison with that of the collapsing fluid. This makes 
the gradient of the scaling law converge to that without the cosmological constant when 
(5 — 5 C ) is sufficiently small. However, the cosmological constant still makes a relevant 
(and growing) contribution on larger scales where the fluid density is lower. 



5. Conclusions 



Following on after our previous work investigating primordial black hole formation 
during the radiative era of the early universe, we have here invested further the critical 
nature of the collapse in this context, focusing on the long-lived intermediate state, 
which appears in collapses of perturbations close to the critical limit, and examining 
the extent to which this follows a similarity solution. This is a key issue, because 
the presence of a similarity solution is seen as an important characteristic feature of 
the general phenomenon of critical collapse. Also, we have presented results from 
calculations where the equation of state parameter w is allowed to take constant values 
different from the radiative value of 1/3, with the aim gaining further insight into our 
main case of interest by considering it within a broader context. Our calculations have 
been made using a purpose-built Lagrangian AMR code, starting with initial supra- 
horizon scale perturbations of a type which could have come from inflation and then 
following self-consistently both the collapse to form the black hole and the continuing 
expansion of the universe. 

From our calculations, we have found that a similarity solution does indeed arise in 
this context, to good approximation, for a region contained within the past light-cone of 
the forming black hole. The self-similarity is not exact, however, and this is interpreted 
as being due to the presence within the light-cone of some outer matter still coupled 
to the expanding universe. We have studied values of w different from 1/3 (up to 
0.6) and found that, for cosmological-type perturbations, scaling-law behaviour persists 
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down to the smallest masses that we are able to follow in each case, with the values 
obtained for the exponent 7 being in almost perfect agreement with those obtained 
previously by Neilsen and Choptuik [15J from their imploding-shells calculation. The 
critical threshold amplitude S c , the intermediate state compactness (2M/R) pea f. and the 
scaling-law exponent 7 all vary with w in a way which is easily understandable in view 
of w being the ratio between the pressure and energy density of the fluid. 
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Appendix: Black-hole and cosmological horizons 

In this appendix we discuss how the concepts of black-hole and cosmological horizons 
emerge from our treatment of light cones. As in the rest of the paper, we assume 
spherical symmetry and that the medium is a perfect fluid. However, the discussion is 
general in the sense that it is independent of the equation of state used and we make no 
assumptions of homogeneity (with reference to the cosmological case), or of asymptotic 
flatness and the presence of vacuum exteriors (with reference to the black holes). We 
think that it is interesting to see how well-known results emerge in this approach. 

We present the argument with cosmic-time slicing, because it is most transparent 
when done in this way, but the final conclusions are independent of the slicing and can 
equally well be demonstrated with the observer time. 

First, we consider the conditions for radial null rays (analogous to the discussion 
in Section 3 for observer time). Recalling that the cosmic-time metric is 



where the plus corresponds to an outgoing null ray while the minus is for an ingoing 
one. The general expression for changes in R along a radial worldline is 



ds 2 = -a 2 dt 2 + b 2 dr 2 + R 2 (d9 2 + sin 2 9d<p 2 ) , 
along radial null rays (for which ds = d6 = dip = 0), we have 



(Al) 



dr = ±-dt , 
b 



(A.2) 



dR 



OR 
~dt 



dt + 



OR 

dr 



dr , 



(A.3) 
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and, using the definitions for U and T given in expressions (|2])-([5]), this can be written 



as 



dR = aUdt + bTdr . 



(A.4) 



Then, inserting condition (|A.2|) into (A.4), we get the expression for how R changes 

(A.5) 



with time along a radial null ray: 
dR , TT _,. 

- = « (C /±r), 

where the plus is for an outgoing ray while the minus is for an ingoing one. Note that we 
are working within the co-moving frame of the fluid and so "outgoing" and "incoming" 
are defined with respect to that frame. 

Next, we recall where the constraint equation ( 13 ) comes from, giving our expression 

(A.6) 



for T. The Gq and G\ components of the Einstein equation give 



4:7rR 2 eR, r = -(R + RU 2 - RT 



and 



AixR 2 apU 



(R + RU 2 -RT 2 ), t , 



(A.7) 



and one can see the same expression appearing in the parentheses on the right hand 

side in each case. It is convenient then to make the definition 
1 



M 



(R + RU — RT ) . 



Integrating equation (|A.6|) gives 
M = 



(Al 



4irR 2 edR, (A.9) 

(corresponding to the interpretation of M as the mass contained within radius iQ, 
while (A.7) gives 

D t M = —AnR 2 pU , (A. 10) 

(corresponding to the change of energy resulting from work done against pressure during 
an adiabatic expansion or contraction). Rearranging the terms in (A. 8) then gives the 

2M 



constraint equation (13), 

r 2 



1 + u 2 



R 



(A.ll) 



Returning now to equation (A.5), the limiting surface at which an outgoing radial 

ren by 

(A.12) 



light ray (ie outgoing with respect to the matter) cannot move to larger R, is given by 
'dR\ 
~d~t ) 



a(U + T) = 0. 



implying 



out 



-u. 



(A.13) 



| Note that there is an effective cancellation here between T factors in both numerator and denominator 
(see Misner [13]). 
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This corresponds to the apparent horizon of a black hole. Similarly, the limiting surface 
at which an incoming radial light ray (with respect to the matter) cannot move to 
smaller R is given by 

Ar a(u ~ r)=0 ' (A - 14) 

implying 

T = U. (A.15) 



This corresponds to the cosmological horizon. Conditions (A. 13) and (A.15) are 
different, of course, but both correspond to 

T 2 = U\ (A.16) 



and then, using (A.ll), to 

R = 2M. (A. 17) 

As noted above, this familiar result (equivalent to that in the Schwarzschild solution) is 
independent of the slicing used and does not require assumptions of homogeneity (with 
reference to the cosmological case), or of asymptotic flatness and the presence of vacuum 
exteriors (with reference to black holes). 
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